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Background: Nuclear pasta phases, present in the inner crust of neutron stars, are associated 
with nucleonic matter at sub-saturation densities arranged in regular shapes. Those complex phases, 
residing in a layer which is approximately 100 m thick, impact many features of neutron stars. 
Theoretical quantum-mechanical simulations of nuclear pasta are usually carried out in finite 3D 
boxes assuming periodic boundary conditions (PBC). The resulting solutions are affected by spurious 
finite-size effects. 

Purpose: In order to remove spurious finite-size effects, it is convenient to employ twist-averaged 
boundary conditions (TABC) used in condensed matter, nuclear matter, and lattice QCD applica¬ 
tions. In this work, we study the effectiveness of TABC in the context of pasta phases simulations 
within nuclear density functional theory. 

Methods: We perform Skyrme-Hartree-Fock calculations in three dimensions by implementing 
Bloch boundary conditions. The TABC averages are obtained by means of Gauss-Legendre integra¬ 
tion over twist angles. 

Results: We benchmark the TABC for a free nucleonic gas and apply it to simple cases such as 
the rod and slab phases, as well as to more elaborate P-surface and gyroldal phases. 

Conclusions: We demonstrate that by applying TABC reliable results can be obtained from 
calculations performed in relatively small volumes. By studying various contributions to the total 
energy, we gain insights into pasta phases in mid-density range. 

PACS numbers: 26.60.Gj,21.60.Jz,02.60.Lj,71.10.Ca 


I. INTRODUCTION 

Nuclear matter as present on earth in the center of 
atoms is almost isotropic with a central density of psat ~ 
0.16 fm“^, the nuclear saturation density. This changes 
drastically in astrophysical environments such as neutron 
stars or core-collapse supernova. In particular, in the 
inner crust of neutron stars at sub-saturation densities 
O.lpsat < P < 0-8/Osati nucleonic matter is expected to 
form complex structures commonly referred to as “pasta” 
phases mm- Because of low proton fractions and macro¬ 
scopic dimensions, pasta phases represent a unique envi¬ 
ronment, which is not present on earth and cannot be 
recreated in the laboratory. 

The occurrence and topology of pasta structures can 
have multiple influences on the properties of a neutron 
star and the evolution of a supernova. Not only it affects 
the neutrino transport but also the neutron star 

cooling [7] and r-mode instabilities in rotating neutron 
stars [SI [3]. 

Nuclear pasta simulations can be divided into two cate¬ 
gories. The first group represents semi-classical methods; 
it includes approaches such as the classical molecular dy¬ 
namics m, Thomas-Fermi approach IIIHI3], and quan¬ 
tum molecular dynamics (QMD) |14H18| . The second 
family includes quantum-mechanical simulations based 
on nuclear density functional theory [1^ , such as Hartree- 
Fock (HF) calculations (with or without BCS pairing) 


|20II27| . While current advanced QMD calculations can 
be performed with hundreds of thousands of nucleons in 
3D boxes greater than 100 fm in size, self-consistent HF 
calculations are still limited to ^2,000 fermions in box- 
sizes of the order of 20 fm. 

Since the large-scale QMD simulations yield periodic 
pasta phases, it is reasonable to assume periodic bound¬ 
ary conditions in HF calculations. In most applications, 
the HF wave functions are also constrained to be periodic. 
This can lead to severe restrictions: Firstly, the solution 
has no chance to develop disorder on a scale larger than 
the box length |2S|- This can only be remedied by us¬ 
ing much larger volumes. Secondly, the periodic ansatz 
for the wave functions is fairly restrictive, because the 
most general solutions for a periodic potential are Bloch 
waves, which differ by a phase when moving to a neigh¬ 
boring box. If strictly periodic wave functions are con¬ 
sidered, spurious finite-size (or shell) effects appear due 
to the quantization of waves in the box. A method to 
remove the spurious finite-volume corrections is to em¬ 
ploy twist-averaged boundary conditions, where the ob¬ 
servables are obtained by averaging over different Bloch 
phases (or twist angles) [29l - l8l] . A method based on 
TABC (sometimes referred to as ‘integration over bound¬ 
ary conditions’) has been applied to circumvent the need 
for large-volume boxes in the context of many-electron 
systems |29II31| . nucleonic matter [32], and lattice QCD 
|331 IMj . It was also used to describe the crust in neu- 
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FIG. 1. (Color online) Nuclear pasta shapes analyzed in this 
work: (a) rod with L = 16 fm, A = 150; (b) slab with L = 
16 fm, A — 294; (c) P-surface L — 22.03 fm, A = 762; and (d) 
gyroid L = 26.01 fm, A = 1254 in a cubic lattice. Dark red 
corresponds to p = 0.14fm“^ and green to p = 0.04 fm“^ 

Iron stars within mean-field models |351155] and Quan¬ 
tum Monte Carlo approach m- The aim of this work is 
to apply TABC to self-consistent HF pasta calculations 
to assess the final-volume effect on previous results ob¬ 
tained with PBC |25fl57] and estimate the realistic box 
sizes with TABC for future applications. 

The manuscript has been organized as follows. In 
Sec. |TTj we describe our implementation of nuclear den¬ 
sity functional theory and twist-averaged boundary con¬ 
ditions. The results of benchmarking calculations for free 
Fermi gas are given in Sec. |III| together with results for 
the rod and slab phases, and triply periodic minimal sur¬ 
face (TPMS) shapes: P-surface and gyroid (see Fig. [^. 
We summarize our results and present the outlook for 
the future in Sec. CYl 


II. METHODS 
A. Skyrme-Hartree-Fock 

To calculate the pasta structures, we use the Skyrme- 
HF method, solving the Schrodinger equation for the 
many-body system in a single Slater determinant ap¬ 
proach. For this purpose we utilize the 3D code Sky3D 
|38| . which solves the self-consistent HF equations with 
a damped gradient iteration method on an equidistant 
grid with no symmetry restrictions. The code, extended 
to the case of TABC, uses the finite Fourier transform 
method for spatial derivatives. We have tested that grid 


spacings between 0.9 fm and 1.05 fm yield stable results 
and this is what is used in this work. 

To calculate the mean field, we use the Skyrme energy 
density functional m-- 

fsk = ^ (cl^{po)p'^ + C^^ptApt 

T=0,1 

+ C^pttt + C^^pt^Jt) , ( 1 ) 

written in terms of the local isoscalar (T = 0) and isovec¬ 
tor (T = 1) densities and currents px (particle density), 
tt (kinetic density), and Jt (spin-orbit current). The 
coupling constants of the functional correspond to the 
Skyrme parametrization SLy6 |39| . as in our earlier work 
|25H27| . The energy density functional is supplemented 
by the kinetic energy and Coulomb terms. 

A uniform electron background is added to ensure 
charge neutrality of the system. Electron screening is 
not included, as its influence should be very small for the 
box lengths used SO]. We consider nuclear matter with 
a proton fraction of Ap = 1/3 to facilitate comparison 
with earlier calculations. 

The computational cost of calculations strongly de¬ 
pends on the number of nucleons and grid size. The 
main time of the calculation is consumed by the diago- 
nalization of the HF Hamiltonian and the wave function 
orthonormalization. Thus, constraining the total aver¬ 
age density is very time-consuming and small systems 
are preferable, at least for the purpose of benchmarking. 

B. Twisted Average Boundary Conditions 

According to the Floquet-Bloch theorem, the single¬ 
particle (s.p.) wave functions for a particle in a periodic 
potential are given by: 

ipaq{r) = Uaq{r)e"'^'^ , ( 2 ) 

where a is a discrete label of the wave function, q is 
the wave vector that determines the boundary condition, 
and Uaq{r) is a periodic function of r. The wave vector 
can be replaced by three twist angles di = qTi {i = 
x, y, z), where Ti = LiSi are the lattice vectors with the 
unit vectors Cp The Bloch waves corresponding to PBC 
represent a particular case of q = 6 = 0. 

The general Bloch boundary conditions can be written 
as: 

tpae{r + Ti) = e^^'tpaeir). (3) 

The s.p.wave functions ipaei'i') defining the HF densities 
and fields are eigenstates of of the HF Hamiltonian hg 
corresponding to the boundary condition (|^: 

hgtpaeir) = €ag'ipae{r). ( 4 ) 

In the TABC method, the expectation value of an ob¬ 
servable O is obtained by averaging over the twist angles: 

( 0 ) = / ( 5 ) 
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FIG. 2. (Color online) Relative finite-size correction for the 
kinetic energy of free symmetric nucleonic matter versus the 
total nucleon number A at fe_F ~ 1.32 fm“^ computed with 
PBS and TABS. The TABS averaging was carried out with 
2, 4, and 10 GL quadrature points in each direction. 

where life is the HF product wave function. In Eq. ([^, 
the angles 9i change between zero (PBC) and tt (anti- 
PBC), as the time-reversal symmetry is assumed |81| . 
In practical calculations discussed in this study, the in¬ 
tegrals over the Bloch phases are computed by means 
of Gauss-Legendre (GL) quadrature with AfGL=4 points, 
unless stated otherwise. 


III. RESULTS 
A. Free particle gas 

Following Ref. m, to test our TABG implementation, 
we consider a gas of non-interacting neutrons and pro¬ 
tons with Xp — 0.5. The exact solution for the average 
kinetic energy can be easily derived from the Fermi gas 
model. Here, we take a system with a total density of 
p = 0.15625 fm“^ {kp « 1.32fm“^) and an average ki¬ 
netic energy of E^in.oo ~ 21.7786 MeV. For plane waves 
in a cubic box L^, the resulting wave numbers are quan¬ 
tized: 

ki,n= -= 0, ±1,±2, ...,±nmax- (6) 

By averaging over 6, we recover a continuous spectrum 
of fe; hence, the infinite-volume limit should be reached 
much faster. 

The results are shown in Fig. The finite-volume 
corrections for different system sizes are plotted for PBG 
and TABS with 2-, 4-, and 10-point GL averaging. The 
results show that the finite-size fluctuations are reduced 
drastically, even with a very small number of GL inte¬ 
gration points. With a higher number of GL points, a 
slightly smoother convergence can be reached. In general, 
with the TABG method, one gains an order of magnitude 
in precision in this case |32| . 


The number of calculations quickly grows as 
where d is the dimension of the phase. For majority 
of pasta phases, it can be assumed that permutations 
of 9i give the same result due to symmetry considera¬ 
tions. For d = 3 the number of calculations decreases to 
{2Nq^ + 6Aql + 4Agl)/ 12. For the following pasta cal¬ 
culations we set A^gl = 4. This means, that we perform 
20 calculations for each system with d = 3. 

B. Pasta phases 

In the following, we apply TABG to actual pasta 
phases. To ensure that the calculations converge at de¬ 
sired shapes, we add a guiding potential ipp during the 


first 

200 HF iterations |27|: 


</'R 

= —(po (cos Y + cos Z ), 

(7a) 

4>s 

= -(po (cosZ ), 

(7b) 

4>p 

= (po (cos X + cos Y + cos Z), 

(7c) 

4>g 

= (po (cos A sin V -|- cos V sin Z -|- cos Z 

sin A), (7d) 


where p G {R, S,P,G} for rod, slab, P-surface, and gy- 
roid, respectively, and Xi = 2 ^x 1 !Li. potentials 

for P and G shapes are their first-order nodal approxi¬ 
mations m- The parameter cpo > 0 has to be chosen 
such that the guiding potential is similar to the resulting 
self consistent potential to guarantee a stable and fast 
convergence. 

1. Rod 

The first example concerns the rod phase. The corre¬ 
sponding shape shown in Fig. [^a) is axially symmetric 
around the cc-axis. We calculate one rod in a rectangular 
box with Ly = Lz = 16 fm. Note that the lowest-energy 
rod configuration corresponds to a honeycomb arrange¬ 
ment m, but the single-rod configuration considered 
here serves well as an illustration of TABG. In this exam¬ 
ple, we vary the particle number A and the box length 
Lx simultaneously to maintain the density and the dis¬ 
tance from the rods in the neighboring cells. The results 
of our calculations are shown in Fig. in the range of 
100 < A < 1200. 

In the limit of large particle number (or L^), the en¬ 
ergy per nucleon should be constant, and this limit is 
reached for particle numbers below A = 1200. In the 
PBG variant, the magnitude of finite-size corrections 
manifesting themselves as energy fluctuations in Fig. 
is large; it reaches « 0.04MeV/A. It is gratifying to see 
that the fluctuations are significantly reduced in TABG. 
For A > 200 the magnitude of finite-size effects falls be¬ 
low O.OlMeV/A. For the potential and kinetic energies 
shown in Figs. ia) and [^b) , respectively, the spurious 
fluctuations are larger; here the TABG method reduces 
the finite-size effects below O.l MeV/A, while they are as 
large as 0.4 MeV/A in the PBG variant. 
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FIG. 3. (Color online) The total energy EjA (bottom), ki¬ 
netic energy i?kin/^ (middle), and potential energy VjA (top) 
per particle for the rod phase of Fig. I^a) at p = 0.0358 fm ® 
computed with PBC (solid line) and TABC (dashed line) as 
a function of particle number A (or the box length L^)- 


2. Slab 

A similar test can be done for the slab phase of 
Fig. gb). This shape is translationally invariant along 
X and y directions. In the first set of calculations we 
keep the box length in the z-direction, = 16 fm, and 
vary the box length in the perpendicular directions si¬ 
multaneously with the particle number to maintain the 
constant density p = 0.0715fm“^. At this density, the 
slab is confined to approximately half the volume of the 
box |27| . 

The results of our test calculations are shown in 
Fig.i Here the improvement provided by TABC is even 
more impressive than for the rod shape, because the 2- 
dimensional averaging is more effective. The plateau in 
E/A is reached well below A = 660 in TABC. Improve¬ 
ments for T/A and V/A are also significant: the range 
of fluctuations is reduced from « l.lMeV/A in PBC to 
« 0.056 MeV/A in TABC. 

Another TABC benchmark for the slab phase can be 
obtained by varying the box length while keeping per¬ 
pendicular lengths constant, = Ly = 16 fm. Here, the 
particle number is adjusted to maintain the thickness of 
the slab at « Tz/2. This has already been done in the 
PBC calculations of Ref. [ 57 ] but the magnitude of finite- 
size fluctuations turned out to be so large that the trend 


Lx=Ly (fm) 
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FIG. 4. (Golor online) Similar as in Fig. but for the slab 
phase of Fig. Sb) at p = 0.0715 fm ® as a function of A (or 

Lx = Ly). 


of the total energy with A was impossible to assess. The 
results of TABC are shown in Fig. (solid line). 

The results for the total energy are shown in Fig.j^a). 
A clear minimum at « 20 fm and E/A = —10.89 MeV 
is found. The reason for this minimum can be understood 
by inspecting different contributions to the total energy. 
To this end, the total energy has been decomposed into 
three parts: the volume energy E^, the surface energy 
E^, and the Coulomb energy Eq- The volume and sur- 


face terms are defined as: 



E'^ — Tfkin + ^ F'Sk.Ti 

T=0.1 

(8) 


E^=Y1 ^Sk.T, 

T=0,1 

( 9 ) 

where 



T^V 

'^Sk,T 

= / (Pr (C'^(po)Pt + EL/pPttt) , 

(10a) 

T^S 

'^Sk,T 

= / <Pr (Ct'’PtA.pt + C^-^pt^Jt) 

.(10b) 

As seen in Fig. [^d), for the small slabs the Coulomb 
energy is very low, because the positively charged slabs 
and the negatively charged voids where the electron gas 
dominates are located very close to each other, and this 
results in a cancellation between electrostatic repulsion 
and attraction. For large slabs, however, the Coulomb 
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repulsion dominates. It appears that this results in an 
almost linear behavior of Ec as a function of L. 



12 14 16 18 20 22 24 26 28 

L(fm) 


FIG. 5. (Color online) TABC energies per nucleon versus 
the periodic length L of the slab shape (solid line), P-surface 
(dashed line) and Gyroid (dotted line) at p = 0.0715 fm“®. 
Shown are: the total energy (a); surface energy (b); volume 
energy (c); and the Coulomb energy (d). 

The volume energy shows a decreasing trend with L; it 
flattens out for large box lengths. At the limit of L —)■ cxd 
the volume energy is supposed to reach the value of « 
— 12.4MeV {E/A for nuclear matter with Xp = 1/3 at 
the equilibrium density of 0.147fm“^). As illustrated in 
Fig. gc), E^ is not close to this value at L = 26 fm, 
because the density within the slab is not yet constant 
and surface effects are still important. 

The surface energy per nucleon in Fig. ib) decreases 
as L~^ = A~^ as the surface area 2LxLy is kept con¬ 
stant and the particle number is proportional to (as 
the particle density is fixed). It is seen that the surface 
energy provides an appreciable contribution to the total 
energy even for the largest box lengths probed in our 
calculations. For small slabs, volume and surface terms 
dominate the behavior of E/A, as both contributions de¬ 
crease rapidly with L. For large box lengths, the pattern 
of E/A is dominated by the Coulomb effect. 


3. TPMS 

At the same mid-density range where the slab phases 
appear, other pasta phases are predicted as well. Those 


are triply periodic minimal surface phases, or TPMS. 
For TPMS shapes, TABC works even better than for the 
slabs, as three-dimensional averaging can be performed. 
Here we vary = Ly = simultaneously to maintain a 
cubic box. The specific shapes analyzed in this work are 
the P-surface shape and the gyroid shown in Figs, [^c) 
and Bd). respectively. Of particular relevance is the gy¬ 
roid phase |421 H5] . Minimal surfaces are of considerable 
interest because of their vanishing mean curvature and 
negative Gaussian curvature. They have been found in 
multiple soft-matter systems |441 - H5] . The nodal approx¬ 
imation of the P-surface and the gyroid are shown in Eqs. 
(7c) and (7dl setting the right-hand side to zero. These 
structures were predicted to appear in nuclear matter in 
time-dependent HF simulations m- 


Furthermore, the P-surface, the gyroid surface, and 
also the D-Surface (diamond; not discussed here), are 
closely related via the Bonnet transformation [13] . Keep¬ 
ing the surface isometric, one can derive the relationship 
between the unit cell lattice parameters: Lp/La = 0.812. 
The ratio of the surface area for the same lattice param¬ 
eters is Ap/Aq = 0.758; this can be compared to the 
slab in the same cubic box: Ap/As = 1.172 jSni- For 
both P- and G-structures, the volume occupied by nu¬ 
clear matter is half of the total volume. Moreover, P, D, 
and G are optimal structures with minimal variations of 
the Gaussian curvature, and G is the optimal structure 
with minimal variations of the structure width 1 501452| . 

The TABG results for the P-surface and the gyroid 
are shown in Fig. At low values of L, gyroids and P- 
surfaces are not stable at certain ranges of twist angles. 
For that reason, we present results only for L > 19 fm. 
The largest calculations for the gyroids in the 29 fm box 
correspond to 1734 particles. This is close to the limit of 
our HF solver Sky3D. 


The minimum of the total energy per nucleon is 
reached at a box length of 21.65 fm for the P-surface and 
at around 28 fm for the gyroid, see Fig. [^a). The surface 
energies in Fig. ib) reflect the ratios of the surface areas 
for gyroid, P-surface, and slab discussed above. Similar 
to the slab case discussed earlier, the surface energy per 
particle for P- and G-shapes decreases as L~^, because 
the particle number is proportional to while the sur¬ 
face area grows as . 

Gompared to the slab case, the Goulomb energy is 
lower for the TPMS at a given box size, because of less 
compact distributions. While for the P-surface this is a 
minor effect, the Goulomb energy for the gyroid is re¬ 
duced by a factor of ~2. This is reflected in the behavior 
of the total energy in Fig.j^a): the minima of TPMS are 
shifted to higher values of L. However, the electrostatic 
repulsion is better compensated by the nuclear energy for 
TPMS; hence, the energy minima for P- and G-shapes are 
not as deep as in the slab case. It should be noted that 
the TPMS minima, especially for the gyroid, are close to 
the energy minimum of the slab phase. 
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IV. CONCLUSIONS 

We implemented the TABC approach into the 3D HF 
framework used to simulate pasta phases in the neu¬ 
tron star crust. We demonstrated that by averaging over 
Bloch boundary conditions one is able to significantly re¬ 
duce the magnitude of spurious finite-volume effects for 
Bloch-Wigner cells containing hundreds of nucleons. 

Practical calculations were carried out for asymmetric 
matter with Xp = 0.3. We first benchmarked TABC for 
the nucleonic gas. The results turned out to be weakly 
sensitive to the number of integration points. We found 
that taking Agl = 4 Gauss-Legendre points yields stable 
results. The TABC method was then tested for the rod 
and slab phases, simultaneously varying the box lengths 
and the number of particles to keep the average density 
constant. For the rod phase, the finite-volume error was 
reduced by a factor of more than three, down to AE = 
0.02 MeV/A. For the two-dimensional slab shapes the 
improvements are even more significant. 

By eliminating spurious finite-size fluctuations through 
TABC, we were able to inspect individual energy contri¬ 
butions from various terms of the energy density func¬ 
tional. This was done by varying the physical length of 
the box. We showed that the energy variation primarily 
comes from the Coulomb and the surface terms. We have 


also demonstrated that, due to its unique geometry, the 
gyroid geometry minimizes the Coulomb energy drasti¬ 
cally for a given box length as compared to the slab and 
the P-surface phases. 

Future applications will include the TABC extension 
of the adaptive multi-resolution 3D Hartree-Fock solver 
[55] and Hartree-Fock-Bogoliubov TABC applications to 
superfiuid pasta phases and complex nucleonic topolo¬ 
gies as in fission. The lesson learned form the exer¬ 
cise presented in this study is that high-fidelity results 
for pasta phases can be obtained by considering finite- 
volume boxes containing up to several thousand parti¬ 
cles. Moreover, TABC calculations are very well suited 
to parallel computing as HF computations at different 
twist angles can be carried out independently. 
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